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ABSTRACT 


The  numerical  method  of  characteristics  with  constant  time  grids 
is  applied  to  the  plane  blast  wave  problem.  The  calculation  begins  at 
a  constant  time  line  bounded  by  the  shock  front  and  a  fixed  back  boundary 
which  lies  halfway  between  the  shock  front  and  the  time  axis.  Along 
the  initial  time  line,  all  three  flow  variables  are  prescribed  and  on  the 
back  boundary,  the  local  flow  velocity  is  prescribed.  The  prescribed 
values  are  calculated  from  the  similarity  solution.  For  the  shock  front, 
the  strong  shock  equations  are  applied. 

When  compared  to  the  exact  similarity  solution,  the  results  of  this 
method  are  found  to  be  accurate  to  within  .5%  for  all  variables  throughout 
the  region  of  calculation,  after  a  pressure  drop  of  98%,  Both  h2~type  and 
h-type  extrapolations  of  calculated  results  improve  the  accuracy;  however, 
h-type  extrapolation  is  more  favorable. 

The  effects  on  the  domain  of  dependence  and  range  of  influence  are 
studied  briefly. 
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NOMENCLATURE 


c  ■  sound  speed 

E  *  parameter  proportional  to  the  total  energy  within  the  wave 

h  »  parameter  proportional  to  the  mesh  size 

k  «  xt/xs  position  ratio  of  back  boundary 
p  ■  pressure 

t  ■  time  variable 

u  =  particle  velocity 

U  «  shock  wave  velocity 

x  *  space  variable 

y  «  specific  heat  ratio 

p  »  density 

Pj  *  constant  density  outside  wave  zone 
Subscripts 

s  *  at  the  shock  front 

b  »  at  the  back  boundary 
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INTRODUCTION 

This  report  is  a  continuation  of  the  report  "Solution  of  Blast 
Wave'  by  the  Method  of  Characteristics."1  To  avoid  unnecessary  repetition, 
the  background  in  the  choice  of  the  particular  problem,  the  governing 
equations,  and  the  similarity  solution  of  the  blast  waves  are  not 
included  in  this  report. 

In  applying  the  numerical  method  of  characteristics,  two  techniques 
are  available;  namely,  the  standard  technique  of  following  the  natural 
characteristic  grid  and  the  technique  of  using  constant  time  increments. 

In  the  standard  technique,  as  developed  in  Ref.  1,  the  points  of  inter¬ 
section  of  the  two  families  of  characteristics  do  not  occur  at  equal 
intervals  in  either  the  space  or  time  coordinates.  When  the  spatial 
distribution  of  the  flow  variables  at  a  later  instant  is  required,  we 
have  to  perform  a  two-dimensional  interpolation,  which  is  not  always 
desirable. 

In  the  present  report,  the  constant  time  technique  will  be  studied. 
This  technique,  originally  proposed  by  Hartree2,  overcomes  the  difficulty 
of  uneven  distribution  of  grid  points  by  choosing  grid  points  on  pre¬ 
determined  constant  time  lines.  This  has  the  advantage  that  the  required 
interpolation  is  always  one-dimensional.  Also,  the  present  method  gives 
results  along  particle  paths  in  successive  steps.  Therefore,  it  is 
advantageous  in  dealing  with  problems  where  interfaces  are  involved. 

Calculations  are  performed  for  the  plane  blast  wave  problem  using 
the  constant  time  method,  and  the  results  are  compared  with  the  exact 
similarity  solution.3  Within  the  region  of  calculation,  indicated  in 
Fig.  3a,  the  results  were  found  to  be  accurate  to  within  .5%,  The  maximum 
pressure  at  the  last  time  line  is  about  2%  of  its  original  value. 

1 


It  sv 


<iy~ 


"•f,  r,~i  »Ti  - 


Three  sets  of  calculations,  each  with  a  different  mesh  size, 
were  made  for  the  same  problem.  The  results  were  extrapolated  with 
h-type  and  h2-type  extrapolation  formulas.  Both  extrapolation  results 
improve  the  accuracy  of  the  results.  The  h-type  extrapolation  gives 
a  better  result. 

It  is  felt  that  the  constant  time  method  will  introduce  error  due 
to  broadening  of  the  domain  of  dependence.  Brief  investigation  indi¬ 
cated  that  this  is  true;  however,  the  error  is  relatively  small. 


FORMULATION  OF  THE  PROBLEM  FOR  NUMERICAL  CALCULATION 

The  one-dimensional  blast  wave  solution  is  one  of  the  few  exact 
solutions  for  the  non-isentropic  (entropy  changes  from  particle  to 
particle),  unsteady  flow.  This  type  of  flow  will  be  calculated  by  the 
numerical  method  (constant  time)  and  the  results  compared  with  the 
closed  form  similarity  solution. 

The  region  in  the  physical  plane  to  be  calculated  is  shown  in  Fig.  1. 

On  the  bottom  it  is  bounded  by  the  constant  time  initial  line.  Along 
this  line,  t  «  tQ,  values  of  all  the  three  dependent  variables  u,  c,  and 
p  are  specified.  In  applying  the  constant  time  method  it  is  desirable 
to  avoid  the  origin  (x  ■  0,  and  t  «  0)  and  the  time  axis.  The  origin 
is  a  singular  point  with  unbounded  values  of  pressure,  particle  velocity, 
sound  speed,  and  along  the  time  axis,  the  sound  speed  is  infinite.  For 
reasons  to  be  given  subsequently  the  time  increments,  At,  between  successive 
constant  time  lines  are  so  chosen  that  they  are  inversely  proportional  to 
the  local  sound  velocity;  in  regions  with  extremely  high  sound  speed  At 
approaches  zero,  making  the  method  impractical.  Hence,  a  back  boundary 
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with  the  equation  xfe  ■  (l/2)xs  is  introduced  on  the  left.  Here,  x^ 
and  xg  are  the  corresponding  x  coordinates  of  the  back  boundary  and 
shock  front,  respectively.  According  to  the  terminology  used  by 
Courant  and  Friedrichs,4  the  back  boundary  is  a  ’’time  like”  arc,  since 
two  of  the  three  characteristics,  II  and  III,  reach  this  curve  from 
the  interior  of  the  region  R.  The  correct  boundary  condition  requires 
that  one  of  the  three  dependent  variables  be  specified  on  the  back 
boundary.  In  this  problem,  the  particle  velocity  u  was  calculated 
from  the  exact  solution  and  used  as  the  one  required  boundary  condition. 
The  right  boundary  of  the  region  of  numerical  solution  is  formed  by 

the  shock  front.  At  a  point  on  the  shock  front,  there  are  four 

unknown  quantities  to  be  determined,  i.e.,  u,  p,  cD  and  U.  The  three 
strong  shock  conditions  together  with  the  equation  of  the  ^character¬ 
istic  serve  to  determine  the  four  unknowns.  The  path  of  the  shock 
front  may  be  established  by  integrating  the  shock  velocity  U. 

NUMERICAL  CALCULATIONS 

Along  the  initial  time  line,  t  *  tQ,  a  number  of  equally  spaced 
points  are  chosen,  and  the  properties,  u,  c,  and  p  at  these  points  are 
calculated  from  the  exact  solution.  Then  a  new  time  line  can  be  drawn 
with  a  time  increment  At.  The  magnitude  of  the  time  increment  is 
limited  by  the  condition  At  <  Ax/c,  where  Ax  is  the  spatial  interval 
between  any  two  neighboring  points  on  the  initial  time  line  and  c  is 
the  corresponding  local  velocity  of  sound.  With  this  criterion,  the  I 

and  II  characteristics  from  a  point  on  the  t  *  t  +  At  line  will  not 

encompass  more  than  one  mesh  point  on  the  t  »  t  line.  Draw  III 
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characteristics  (path  iine)  from  each  of  the  initial  points »  The 
intersections  of  these  path  lines  and  the  new  time  line  establish  a 
set  of  new  points  about  which  the  flow  properties  are  to  be  evaluated. 

To  accomplish  this,  the  state  and  physical  characteristic  equations 
have  to  be  solved.  Since  the  equations  for  physical  characteristics 
involve  the  state  variables  u  and  c  and  the  state  characteristic 
equations  contain  the  physical  variable  t,  neither  of  them  can  be  solved 
independently.  Numerical  iterative  schemes  are  devised  to  obtain 
approximate  solutions  and  the  detailed  iterative  procedures  are 
described  in  Appendix  A.  Referring  to  Fig.  A. 2,  the  iteration  procedure 
involves  establishing  the  Ill-characteristic  through  an  initial  point  2, 
locating  the  intersection,  i.e.,  point  4,  of  the  particle  path  and  the 
constant  time  line  t  *  t  ♦  At,  tracing  back  the  I  and  II-characteristics 
from  point  4  until  they  reach  the  initial  time  line  at  points  A  and  B. 

The  properties  at  A  and  B  are  calculated  by  interpolation  of  values  at 
points  1,  2,  and  3.  The  fact  that  the  entropy  along  a  particle  path 
remains  constant  is  used  in  relating  properties  at  point  2  and  point  4. 

The  form  of  the  interpolation  formulas  used  in  obtaining  properties 
at  points  A  and  B  are  highly  influencial  to  the  accuracy  of  the  constant 
time  method.  Linear  interpolation  would  be  the  simplest  in  computation. 
Unfortunately,  as  shown  in  Fig.  2,  the  pressure  p  near  the  shock  front 
and  the  sound  velocity  close  to  the  back  boundary  are  both  highly 
non-linear  with  respect  to  the  spatial  coordinate  x.  Linear  interpolation 
would  cause  appreciable  error  for  the  values  of  these  properties. 
Furthermore,  the  errors  tend  to  accumulate  with  time.  Errors  due  to 
linear  interpolation  which  might  be  tolerable  in  usual  cases  are 
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objectionable  in  the  present  problem.  For  this  reason  quadratic  interpolation 
formulas  were  used  in  interpolating  p  and  c  while  linear  interpolation  was 
used  for  u.  Local  velocity  u  is  essentially  linear  with  respect 
r,o  x  throughout  the  region  of  solution  R. 

At  the  back  boundary  and  the  shock  front  modified  iteration 
procedures  are  adopted.  Since  the  shock  front  propagates  outward 
engulfing  new  particles  continuously  and  in  the  mean  time  flow  particles 
leave  the  region  of  solution  across  the  back  boundary,  sub-routines  are 
provided  to  terminate  a  particle  path  which  is  about  to  leave  and  to 
initiate  a  new  particle  path  near  the  shock  front  when  the  distance 
between  the  shock  front  and  the  closest  data  point  equal  or  larger  than 
the  local  Ax.  Therefore,  the  number  of  data  points  between  successive 
time  lines  varies.  In  the  present  case,  the  total  number  of  data  points 
along  a  constant  time  line  increases  slightly  with  time.  This  means 
there  are  less  particle  paths  passing  through  the  back  boundary  than  the 
new  particle  paths  initiated  at  the  shock  front.  The  rate  of  increase 
of  data  points  is  slow;  which  is  not  enough  to  cause  concern  about 
consuming  too  much  computation  time* 

The  numerical  procedure  was  programmed  in  Fortran  II  language  and 
the  computation  was  done  on  an  IBM  7040  computer.  The  computation  time 
for  the  case  of  51  initial  data  points,  h  *  1/2,  was  approximately  half 
an  hour,  which  gave  a  pressure  drop  of  98%  on  the  shock  front.  The 
accuracy  of  this  calculation  is  shown  in  Fig.  3.  Figure  3a  shows  the 
physical  plane,  indicating  the  position  of  the  shock  front,  the  back 
boundary,  and  the  last  line  of  numerical  calculation.  Figure  3b  shows 
the  relative  errors  of  p,  u,  and  c  at  points  on  the  last  line  of  the 
calculation  as  compared  to  the  exact  solution.  The  relative  errors  of 
t,  p,  u,  and  c  along  the  shock  front  are  shown  in  Fig.  3c.  The  maximum 
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error  for  all  variables  throughout  the  entire  region  of  calculation 
is  less  than  .51. 

EXTRAPOLATION  OF  NUMERICAL  RESULTS 

In  order  to  study  what  type  of  extrapolation,  h-type  or  h2-type 
(Ref.  1),  may  improve  the  numerical  results,  the  same  plane  blast  wave 
problem  was  calculated  three  times  using  three  different  numbers  of 
initial  data  points.  Identical  computations  were  performed  with  26,  51 , 
and  101  initial  points  which  correspond  to  mesh  sizes  h  *  1,  1/2,  and 
1/4,  respectively.  The  relative  errors  in  p,  u,  c,  and  t  along  the 
shock  front  for  the  three  cases  are  plotted  in  Fig,  4,  The  curve 
labeled  h  ■  1/2  represents  the  error  of  the  indicated  properties  with 
mesh  size  h  «  1/2 „  The  errors  decrease  as  the  mesh  size  decreases  and 
the  general  shape  of  these  error  curves  resemble  each  other  within  their 
own  group.  The  highest  error  appears  in  p;  while  the  least  in  u  and  c. 
The  extrapolated  values  of  p  along  the  shock  front  were  compared  with 
the  exact  solutions.  Results  are  presented  in  Fig,  5  and  Fig,  60  In 
these  figures  the  original  errors  were  also  plotted  for  easy  comparison. 
The  curve  labeled  h2(l,  1/2)  represents  the  error  in  the  extrapolated 
value  of  p  calculated  from  the  h2-type  formula,  by  using  the  pressures 
from  corresponding  calculations  with  h  =  1  and  h  *  1/2;  h(l,  1/2,  1/4) 
represents  h-type  extrapolation  of  pressures  corresponding  to  h  »  1, 

1/2,  and  1/4,  etc. 


These  figures  indicate  that: 

14  The  extrapolated  results  from  any  two  sets  of  calculations  are 
always  more  accurate.  This  is  true  for  both  h-type  and  h2-type 
extrapolation, 

2.  Results  obtained  from  the  h-type  extrapolation  are  more  accurate 
than  those  from  the  h2-type  extrapolation. 

3.  Extrapolation  from  three  sets  of  calculated  results  are  slightly 
better  than  those  from  two  sets  of  calculations. 

DOMAIN  OF  DEPENDENCE  AND  RANGE  OF  INFLUENCE 

It  is  known  from  theory  of  partial  differential  equations  of  hyper¬ 
bolic  type  that  the  domain  of  dependence  of  a  point  in  the  physical  plane 
lies  on  the  portion  of  the  initial  line  that  is  bounded  by  the  I  and  II 
characteristics  from  the  point.  For  instance,  any  change  in  flow 
properties  outside  the  domain  AB  in  Fig.  A.2,  certainly  will  not  effect 
the  properties  at  point  4  theoretically.  By  using  the  constant  time 
method,  properties  at  points  A  and  B  are  determined  by  interpolating 
properties  at  points  1,  2,  and  3.  Consequently,  any  change  of  values 
of  properties  at  points  1  and  3  will  in  fact  change  the  values  of 
properties  at  points  A  and  B,  and  in  turn  will  change  the  values  of  the 
properties  at  point  4.  Apparently  this  is  an  error  attributed  to  the 
numerical  method.  The  criterion  At  *  R(Ax/c),  where  R  <  1,  used  in 
the  choice  of  At  is  equivalent  to  requiripg  that  points  A  and  B  lie 
inside  of  point  1  and  3,  The  broadening  of  the  domain  of  dependence 
from  AB  to  73  can  also  be  viewed  as  the  widening  of  the  range  of 
influence  of  a  point  on  the  initial  time  line. 
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To  evaluate  the  error  introduced  by  this  broadening  of  the  domain 
of  dependences  several  calculations  were  made  for  the  same  problemc 
with  computations  confined  in  a  small  region  R  in  Fig0  7a.  The 
prescribed  properties  for  the  equally  spaced  points  on  the  initial 
line  were  calculated  from  the  exact  solution.  The  method  of  character** 
istics  was  first  used  in  the  calculation  to  determine  the  II~character~ 
istic  passing  through  point  L0  This  established  a  boundary e  shown  in 
Fig.  7b.  To  the  left  of  this  boundary p  the  flow  properties  should  not 
be  influenced  by  any  changes  in  properties  at  and  to  the  right  of 
point  L.  Then  four  sets  of  calculations  were  performed  using  the 
constant  time  method.  The  first  set  to  be  used  as  a  standard  for  com^ 
parison  is  calculated  from  correct  prescribed  initial  values.  In  the 
other  three  calculations  the  flow  properties  at  point  L  were  changed. 

The  changes  were  made  by  varying  one  of  the  three  properties e  uc  cp 
or  p,  and  keeping  the  other  two  unchanged.  For  example^  in  the 
calculation  involving  a  change  in  p0  its  value  was  reduced  by  5%  at 
point  L.  The  deviation  of  the  pressure  from  the  first  set  of 
calculation  in  the  region  to  the  left  of  LP  were  calculated  and  indicated 
in  Fig.  7b„  Similar  results  for  u  and  c  may  also  be  obtained.  Such 
deviations  are  deemed  as  errors  due  to  the  constant  time  method.  It 
should  be  mentioned  here  that  the  time  increments  in  all  the  four  sets 
of  calculations  were  the  same  but  the  x  coordinates  of  the  data  points 
do  not  coincide  exactly  between  different  sets  of  calculations.  In 
computing  the  relative  errors 9  such  deviations  in  x  coordinates  were 
compensated  through  linear  interpolation.  The  range  of  influence  of 
point  L  is  widened  and  it  is  bounded  on  the  left  by  LQ  which  has  a  slope 
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of  (dx/dt)  *  u  -  (c/R)  rather  than  (dx/dt) jj  >  u  -  c.  Since  R  is 
always  less  than  one,  the  dotted  line  will  therefore  be  situated  to 
the  left  of  the  II-characteristic,  LP,  However,  the  error  is  compara¬ 
tively  small  and  tends  to  diminish  gradually  with  time. 

CONCLUSION 

The  jnumerical  method  of  characteristics  .with  constant  time  grids 
produces  very  accurate  results  when  applied  to  plane  blast  wave  problem, 
provided  the  region  with  extremely  high  sound  velocity  is  excluded.  For 
other  one-dimensional  unsteady  fluid  flow  problems  which  do  not  include 
regions  with  extremely  high  velocity  of  sound,  it  is  reasonable  to 
believe  that  this  numerical  method  will  also  yield  accurate  results. 

Both  h-type  and  h2-type  extrapolation  >f  the  numerical  results 
improve  the  accuracy,  whereas  h-type  extrapolation  gives  a  better  result. 
The  accuracy  of  the  numerical  results  is  sensitive  to  the  interpolation 
formula  employed.  This  method  is  particularly  useful  for  problems 
involving  interfaces  which  need  to  be  traced. 

Although  the  results  of  the- plane  blast  wave  problem  only  are 
reported,  the  Sana  method  with  slight  modification  can  be  used  to  solve 
problems  with  cylindrical  and  spherical  symmetry.  Work  is  being  done  in 
this  area  and  will  be  reported  at  a  later  date. 
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Figure  1  Physical  Plane  Showing  Region  of  Numerical  Solution  for  Blast 
Wave  Problems 


t./sec.  x  10“ 


I 

Figure  3  Plane  Wave  Error  Analysis 
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APPENDIX  A 


Iteration  Procedure 

In  this  appendix-  the  detailed  iteration  procedure  for  the  constant 
time  method  is  described. 

The  basic  step  in  the  solution  of  the  unsteady  one^dimensionai 
fluid  flow  problem  by  the  constant  time  method  involves  the  determination 
of  all  flow  variables  u,  cs  pt  x  at  a  set  of  discrete  points  along  a 
constant  time  line,  when  those  variables  are  known  along  a  previous 
time  line-.  As  it  was  mentioned  before  the  choice  of  at  is  subject  to 
the  condition  at  <  Ax/c  for  all  interior  points .  except  points  on  the 
back  boundary  and  those  on  the  shock  front.-  For  this  particular 
problem  we  are  fortunate  to  be  able  to  determine  the  proper  at  based  on 
properties  near  the  back  boundary  without  searching  through  the  complete 
initial  time  line  More  specifically,  the  criterion  used  (ref;  Fig:  A  1) 
is,- 

At  =  w^ere 

Repeated  application  of  the  above  step,  with  proper  boundary  conditions, 
yields  solution  of  the  problem  For  the  blast  wave  problem,  we  encounter 
tour  basic  types  of  points  They  are; 

1.  Points  on  the  back  boundary 

2.  Interior  points  which  lie  within  the  region  cf  numerical  solution 

3  Points  on  die  shock  tront 

4.  Points  on  added  path  line 
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Basic  unit  operations  used  in  solving  these  four  different  types  cf 
points  are  described  as  follows? 


1,  Points  on  the  Back  Boundary 


Referring  to  Fig,  A„la  all  variables  (u,  c,  p,  x,  t)  at  points  1,  2,  and 
3  are  known,  the  variables  at  point  4  are  to  be  determined.  Points  1 
and  4  are  on  the  back  boundary.  Points  4  and  A  are  on  the  same  III 
characteristics  while  points  B  and  4  are  joined  by  a  common  I I -charset er¬ 
istic,  Point  4  must  satisfy  the  equation  of  the  back  boundary  which  is 

*4=  <)V  <*•» 

where  k  *  constant  (the  position  ratio) ,  and  as  discussed  earlier,  the 
particle  velocity  along  the  back  boundary  is  prescribed  by  the  exact 
solution 


F3<r+0  ( f, 


JL 

3 


(A.  2) 


19 


where  F  *  constant,  (the  velocity  ratio  u^/u^.  Since  At  has  been 
determined,  t 4,  the  time  at  point  4,  is  known  to  be  tQ  ♦  At.  Therefore 
x4  and  u4  can  be  calculated  readily  from  equations  (A.l)  and  (A.2), 
respectively.  The  finite  difference  forms  of  the  physical  character* 
istic  equations  are: 


Along  II 


At 


(A.3) 


Along  III  AA  =  (A.4) 

where  the  barred  quantities  represent  the  average  values,  i.e., 


The  corresponding  state  characteristic  equations  are: 


Along  II 

fV-P. 

a  r-5t 

(A.S) 

Along  III 

ii 

Ofjvf 

*  Fa) 

(A. 6) 

The  properties  at  points  A  and  B  are  interpolated  from  the  known  values 
at  points  1  and  2,  or  2  and  3.  Since  u  is  very  much  linear  with  respect 
to  x  throughout  the  region,  linear  interpolation  formulas  are  used  to 
determine  u  at  points  A  and  B. 


U,  "t~  38  fy**) 
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For  the  other  properties s  c  and  pp  at  points  A  and  Bp  quadratic 
interpolation  formulas  are  applied,, 


A JXt-xJ+Bctez-Xt)  TCl  =  C(*B ) 

(Ac  9) 

Pg  ~  Ap  CXg“Pt)  +  &f(Xg  ^  Pg  —  F(Xg) 

(A* 10) 

Ca  =  At(XA-Xtft3c(XA-x,)-h  C4  =  CCxA) 

(A.  11) 

Pa  =  Ap  (Xa  -Xif*  &?(Xa-'QT  Px  =  F(*a) 

(A. 12) 

where 

A,-‘: 

(Xs-bf  *,-** 


?c  = 


&-*d\  frf* 

PET 

The  ten  equations  (A„ 3)  to  (A. 12)  can  be  solved  for  the  ten  unknowns 
XA®  xbp  uAs  ubp  cAC  cBP  pAP  pBP  c4#  and  p4,  The  solution  is  accomplished 
by  the  following  iteration  scheme,,  To  avoid  possible  confusion  in 
notations  in  the  description  of  the  iteration  procedures p  a  convention 
in  using  the  superscripts  will  be  adopted  throughout  this  appendix  unless 
otherwise  specified,,  Superscript  os  on  any  variables p  such  as  uA„ 
represents  the  initial  estimation  of  the  value  at  point  A„  A  single 


C3~C* 

PET 

frfk  *r* 
r3-r* 

PET 


prime  on  any  variable,  such  as  u^,  designates  the  first  approximation 
at  point  B.  Similarly  a  double  prime  on  any  variable  indicates  the 
second  approximation  of  the  value  of  such  a  variable,  etc. 


First  Iteration 

Rewriting  equations  (A. 3)  to  (A, 12)  in  the  above  mentioned  notation 
and  making  some  approximations. 


(A. 13) 

(A. 14) 

(*a) 

(A. IS) 

(A. 16) 

Cj  =  CM) 

(A. 17) 

P»  =  P(*e) 

(A.  18) 

d  =  C«> 

(A. 19) 

Pa  =  ?M) 

(A. 20) 

P*  =  r('7£c{)('**-'A-i)*Pe 

(A. 21) 

C  '  -  (Jl\  *fr  ' 

C*  ~  \  Pa)  -a 

fA.22) 

The  iteration  procedure  is  initiated  by  estimating  values  for  u°}  c°, 
c°.  u°.  and  p°.  ,  It  is  reasonable  to  assume  that 


Ue  =  iAa  =  U; 


* 


z 
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With  these  estimations,,  x®  and  can  be  calculated  from  (A,  13)  and 
(A.  14)  8  respectively,,  Once  x^  and  x^  are  determined  the  properties 

UA°  UB*  CBS  ^Bc  CAD  an<*  ^A  can  comPut:ec*  ^rom  interpolation 
formulas  (A.15)  to  (A. 20).  Finally,,  from  equations  (A.21)  and  (A0 22) £, 

Pi*  and  c°^  may  be  obtained  independent ly0  The  first  approximation  of 
the  properties  at  point  4  and  at  points  A  and  B  is  now  complete . 

Second  Iteration 

Once  the  first  approximation  is  accomplished,,  the  second  and 
subsequent  approximations  of  the  variable  at  point  4  can  be  established 
from  a  set  of  equations  similar  to  equations  (A, 13)  and  (A. 22).  Repeating 
the  procedure  in  the  first  iteration  with  the  primed  values  replacing 
the  estimated  values  (values  with  superscript  "o")  and  double  primed 
variables  instead  of  single  primed  variables,,  we  can  obtain  the  second 
approximation,,  This  iteration  process  may  be  continued  until  the  desired 
degree  of  accuracy  is  reached.  In  the  present  calculation „  the  iteration 
was  stopped  when  the  pressure  at  point  4  converged  to  within  .0005%  of 
its  value  in  the  previous  iteration.  For  most  points  in  th:is  problem,, 
not  more  than  four  iterations  were  required  to  achieve  the  desired 
accuracy. 

2,  Interior  Points 

» 

+  4 1 


FigmcA.7.  i  hysical  Plane 


Referring  to  (A. 2),  all  variables  are  assumed  known  at  points  1,  2,  and 
3,  and  At  is  predetermined.  The  variables  at  point  4  are  to  be 
determined.  Point  A  and  point  B  lie  on  the  characteristics  along  I 
and  II  issued  from  point  4,  respectively,  and  points  4  and  2  are  on  the 
same  path  line  (Ill-characteristic).  Since  At  is  chosen  such  that 
At  <  Ax/c,  therefore  for  most  interior  points,  point  A  lies  between 
point  1  and  2,  and  point  B  between  points  2  and  3.  For  interior  points 
next  to  the  back  boundary  or  shock  front  some  exceptional  cases  occur. 
These  special  cases  are  covered  later  in  this  section. 


The  finite  difference  forms  of  the  physical  characteristic 
equations  are: 


Along  I 

X(.-Xa  -  r 

At  =  ^A4  ■** 

(A. 23) 

Along  II 

At  - 

(A. 24) 

Along  III 

M~Xt  77 

At  -  Hff 

(A. 25) 

The  corresponding  state  characteristics  are: 

Along  I 

P*-f*  /  Pm 

-~r  yi 

(A. 26) 

Along  II 

fr-fb  f  Pa 4 

(A. 27) 

Along  III 

ii 

S’ I*? 

(A, 28) 

The  following  interpolation  formulas  for  properties  at  points  A  and  B 
are  the  same  as  defined  before, 


“a*  Va(*a) 

(A. 29) 

Vw 

p? 

II 

■0 

75 

(A.  30) 

—  C  (X$) 

(A, 31) 
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(A«32) 


r,=  F«g) 

Crf  =  C(*t)  (A.  33) 

Pa  -  PM  ca.34) 

From  the  twelve  equations  „  (A.  23)  to  (A034)s  the  twelve  unknowns 0  x4P  u,^ 
p4«  C4B  XA»  uAe  CA0  PA»  XB*>  UB0  CB8  PB  can  be  solvedo  The  approximate 
solution  is  obtained  by  the  following  iteration  procedure. 

First  Iteration 

The  first  approximation  was  begun  by  rewriting  the  system  of 
equations  (A023)  to  (A. 34)  in  a  new  form  which  involved  some  approximation. 
We  have 


<=  *  1- (**&*)** 

(A„35) 

(A„36) 

xb=4-(^-^M 

(A037) 

u 

'■N 

A 

(A„38) 

Mg  ~  Ug(Xe) 

(A.39) 

Cg  =  c  (Xe) 

(Ao40) 

P6'=  ?(*s) 

(A.41) 

4 =C(Xa) 

(A042) 

3>» 

li 

~X) 

£ 

1 

(A.43) 

?J-Pa  s/ 4+ %  \ 

Uj-U/l  ~  ' 

(A044) 

(A045) 

4  = 

(A046) 
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and  x£  is  determined  by  assuming  u4  =*  U20  Now  ard  may  be  calculated 
from  equations  (A. 36)  and  (A. 37)  with  the  following  estimated  values; 
u°  «  Uj}  u°  =  u3 9  c°  *  cls  Cg  =  c39  c°  =  c2.  Substituting  x^  in  equations 
(A, 38)  s  (A.42)p  and  (A.43)„  the  properties  at  point  A0  uj^  cj^  and  p^„ 
can  be  determined,,  and  substituting  Xg  in  equations  (A.39),,  (A, 40),,  and 
(A. 41)  determines  u’y  c°„  and  p^»  Setting  p°  equal  to  p2£)  equations  (A. 44) 
and  (A.4S)  can  be  solved  simultaneously  for  the  properties  pJ4  and  u'4  0 
Finally  c}4  is  computed  from  equation  (A*46) «,  The  first  approximation 
of  the  properties  at  point  4  and  points  A  and  B  are  now  complete 0 

Second  Iteration 

By  repeating  the  above  procedures  in  the  first  iteration  with 
single  primed  values  replacing  the  estimated  values  and  double  primed 
variables  instead  of  single  primed  variables 0  the  second  approximation 
can  be  established.  The  iteration  process  nay  be  continued  until 
the  desired  degree  of  accuracy  is  reached.  In  the  present  calculation,, 
the  iteration  was  stopped  when  the  properties  u4  and  p4  are  both 
converged  to  within  .0005%  of  their  corresponding  values  in  the  previous 
it-  ^’tion. 


Point  4  lies  on  the  left  of  the  back  boundary.  Since  it  lies 
outside  the  region  of  our  solution,,  computation  for  such  points  are  skipped. 


Figure  A. 2a  Physical  Plane 
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back  boundary 


The  ^characteristic  issued  from  point  4  intersecting  the  back 
boundary  at  A  before  reaching  t  *  tQ0  The  general  iteration  procedure 
for  the  interior  point  was  followed  except  modifications  were  made  to 
determine  the  location  of  point  A  and  the  properties  at  A.  Here  the 
local  back  boundary  between  points  1  and  5  was  approximated  by  a  straight 
line 

Xa  ~Xt  _  Xf—X,  (A„47) 

t/  ■  if  ~  a  t 


The  finite-difference  form  of  the  I -characteristic  equation  in  the 
physical  plane  is 

.  4*  jit  (A.48) 

t4~tA  ~  Z  Z 


After  x'4  is  determined  by  equation  (A,35)s  the  first  approximation  of 
the  location  Ap  x^  can  be  solved  from  the  following  version  of  equations 
(A, 47)  and  (A„48)  (i.eop  equations  (A. 49)  and  (A, SO))  with  the  indicated 


estimated  values B  bearing  in  mind  that  point  5  is  a  known  point  by  this  time. 


XA~ X) _ Ks 

tA-t0  ~  At 

X*-Xa  U-A  +  UA  Ca *j-<4  uA  -  U, ' 

tj-tA  ~  Z  +  z,  >  and  tA-C,J 


(A, 49) 


(A,  50) 
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The  properties  at  point  A  were  determined  by  linear  interpolation 
formulas 


The  II- characteristic  intersects  the  shock  front  before  it  reaches 
the  line  t  =  t  t  Points  next  to  the  shock  front  are  always  computed 
after  the  properties  at  the  shock  front  (point  5)  has  been  determined. 
The  actual  procedure  in  establishing  the  properties  at  point  5  will  be 
described  later.  Here  we  consider  x5i>  t5S  u5i)  c5„  and  p5  to  he  known 
values.  The  modification  in  the  iteration'  procedure  for  this  case  is 
similar  to  that  of  special  case  Bc  After  x54  is  determined  from  equation 
(A/35jp  and  assuming  that  the  shock  front  between  points  3  and  5  to  be  a 
straight  line,,  x£  can  be  solved  from  the  equations v 

D 

xB-X3  Xs 

1B-*0  ~ 

v-V  t  u*°+  u‘  fsel 

{  and  J 
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The  linear  interpolation  formulas*  which  are  used  to  determine  the 


properties  (first  approximation)  at  point  B*  are 
U'= 

Cg  ~  Cj-h  ( 

fs  =  + 


3,  Po*nts  on  the  Shock  Front 


shock  front 


From  the  definition  of  shock  wave  velocity 


shock 


we  have 


X4—X3 

At 


(A, 51) 


The  finite*»difference  form  of  the  I»physical  and  state  characteristics  are* 


-  .7  .pr 

a  V  “  *  CH 

?4~Fa _ y' 

**4  ^4 


(A, 52) 
(A. 53) 


The  properties  at  point  4  must  satisfy  the  strong  shock  equations.  Therefore 

1  - 


u 


C+-J 


2U 

f.(r+i ) 

gcEU 

z, 


(A, 54) 
(A055) 


Linear  interpolation  formulas  aTe  used  to  determine  u  and  c  at  point  Ae 
Since  in  this  case  point  A  may  lie  anywhere  between  point  1  and  point  3, 
two  sets  of  formulas  (A,56a)s  (A* 57a)*  and  (A, 56b)*  (A, 57b)  are  provided,, 
If  A  lies  between  points  1  and  2 


4.-  =£'<*>  (A-57a) 


or  if  A  lies  between  points  2  and  3 


4a“  V4(*a) 

CA  sr  Cz  +"(  Xj-Xj  )  (  S 


(A, 56b) 
(A057b) 


Since  p  varies  rapidly  with  respect  to  x  near  the  shock  front „  the 


quadratic  interpolation  formula  is  applied0  From  equation  (AUI2) 


a-  at( r*  =  Pr* )  (A-58) 

From  the  eight  equations 0  (A. 51)  to  (A.58)c  the  eight  unknowns  x^B  u „ 
ciis  P^s  x^b  ^A£l  cAP  and  pA  may  be  determined,  The  solution  is  accomplished 
by  the  following  iteration  procedure. 


First  Iteration 

Equations  used  in  the  first  approximation  ares 


(A,59) 

(A, 60) 

=  U/(%\)  or  Uz(%a) 

(A„61) 

c\ 

*  . 

II 

Pi 

-TN 

s 

o 

H 

'“N 

(A,62) 

II 

& 

(A063) 
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(A, 64) 


(A065) 
(A,  66) 


where  equations  (A. 61)  and  (A. 62)  are  the  appropriate  linear  interpolation 
formulas.  The  first  iteration  is  performed  with  the  following  initial- 
input  values  %  u°  =  u3P  c°  *  c3 p  p°  *  p3 #  u°  »  u2fi  and  c°  *  c2.  The 
first  approximation  of  all  the  eight  variables  can  be  computed  one  by 
one  from  equations  (A. 59)  to  (A.66)0 


Second  Iteration 

The  second  iteration  is  accomplished  by  replacing  the  initial 
input p  c°p  and  p°  (estimated  values) 6  with  the  average  values  of 
the  initial  input  and  the  first  approximation  of  the  corresponding 
variable3  u°  by  u^8  c°  by  c^p  and  replacing  single  primed  variables 
with  double  primed  variables  in  equations  (A.59)  to  (Aa66)p  and  computing 
the  eight  variables  following  the  sequence  of  equations „  Repeating  the 
second  iteration  using  the  average  values  of  the  first  approximation  and 
the  second  approximation  as  inputs  and  substituting  the  double  primed 
variables  with  triple  primed  variables B  we  can  obtain  the  third 
approximation.  The  iteration  procedure  may  continue  until  the  desired 
accuracy  is  reached.  In  this  casep  the  iteration  was  stopped  when  the 
pressure  at  point  4  converged  to  within  ,0002%  of  its  value  in  the  last 
iteration.  The  reason  for  using  the  average  values  of  the  initial  input 
and  the  current  approximation  as  the  new  input  in  the  next  iteration 
rather  than  just  the  current  approximation  of  the  corresponding  variables 
is  because  the  later  scheme  results  diverge  in  subsequent  iterations. 
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Another  iteration  scheme*  has  been  tested  and  compared  with  the  one 
just  described,,  This  scheme  gives  a  convergent  result  with  desired 
accuracy 5  however  it  converges  at  a  slower  pace0 

4*  Points  on  Added  Path  Line 


shock  front 

* 

Figure  A„4  Physical  Plane 

Referring  to  Fig„  At4p  when  Ax,,  >  Ax^  a  new  particle  path  will  be 
initiated  from  the  shock  front  for  the  purpose  of  keeping  the  data  points 
populated  more  or  less  uniformly  throughout  the  region  of  calculation 0 
As  mentioned  in  special  case  C  point  4  is  always  treated  after  point  5 
has  been  established,,  The  present  case  is  exactly  the  same  as  in  special 
case  C  except  point  2  and  point  3  coincide „  Hence „  the  same  iteration 
procedure  can  be  used  to  determine  the  properties  at  point  4„ 

♦Instead  of  computing  p'^  and  u'4  separately  from  (A„64)  and  (A065)  ^  we 
could  solve  the  following  two  equations  simultaneously  for  p54  and  u’. 
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